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Processing Seismic Data 

The present invention relates to processing seismic data and, in particular, to processing 
marine seismic data to reduce the effects of "ghost reflections". The present invention 
can be applied to processing both marine seismic data acquired in a flat sea and marine 
seismic data acquired in a rough sea. 

Figure 1 is a schematic diagram of a marine seismic survey in which seismic energy is 
emitted from a source 1 and detected by a seismic receiver 2 at a depth h below the 
surface 6 of the sea. Energy emitted from the source is reflected by the seabed 3 or by a 
reflector 4 below the seabed 3 and is then detected by the receiver. This path of seismic 
energy is labelled 5 in Figure 1. Information about the geological structure of the 
earth's interior can be derived from the reflected seismic energy incident on the receiver 

The seismic receiver 2 shown in Figure 1 is a streamer, which is a type of seismic 
receiver that is frequently used in marine seismic surveying. A streamer contains a 
plurality of sensors Si, S2,...Sn such as pressure sensors and/or particle velocity sensors 
distributed along its length, which may be some hundreds of metres, and is thus able to 
measure the reflected seismic energy at a number of points simultaneously. A streamer 
may be suspended from one or more floats 8, so that all the receivers of the streamer are 
at the same depth in a flat sea. 

In addition to the desired path 5 of seismic energy shown in Figure 1, other seismic 
energy paths will occur as a result of seismic energy being reflected or scattered from 
the sea surface 6. These paths are known as "ghost reflections". For example, reference 
7 in Figure 1 shows a ghost reflection in which the seismic energy reflected by the 
reflector 4 is not directly incident on the receiver 2, but undergoes a further reflection at 
the sea surface 6 before reaching the receiver. Down-going sea-surface ghost 
reflections are an undesirable source of contamination of seismic data, since they 
obscure the interpretation of the desired up-going reflections from the earth's interior. 
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As is well-known, the ghost signal will produce one or more "notches" in the frequency 
spectrum of seismic energy, with the frequency at which the notches occur depending 
on the depth of the receiver below the sea-surface. 

The ghost signal reflected from the sea surface is delayed relative to the direct signal. 
There are two components to this delay: firstly, there is a 180° phase change upon 
reflection at the sea surface and, secondly, there is a time delay corresponding to the 
additional path length (which for a signal emitted in the vertical direction is 2z, or twice 
the depth of the receiver). The actual vertical far field signal is the sum of the direct 
signal and the ghost signal. The direct signal and the ghost signal will interfere, and this 
causes variations in the amplitude of the far field signal. For some frequencies, the 
interference is destructive and causes a zero amplitude or "notch" in the spectrum. This 
occurs at frequencies where the depth of the receiver is an even number of quarter 
wavelengths: 

./notch = (nc/2z\ n - 0, 1, 2, 3. . . 

(where c is the speed of sound in water, n is an integer giving the harmonic order, and z 
is the depth of the receiver below the sea surface). 

Constructive interference occurs at frequencies exactly intermediate adjacent notch 
frequencies, and this leads to maxima in the amplitude at these frequencies, given by: 

^eak = (2« + \)c/4z, n = 0, 1, 2, 3. . . 

The effect of the interference between the direct signal and the ghost signal can be 
thought of as applying a frequency domain "ghost filter" to the direct signal. The ghost 
filter has the following form: 



£09 = 1+ I r | 2 -2 | r | cos(2z//c) 
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(where r is the reflection coefficient at the sea surface). The amplitude spectrum of the 
ghost filter resembles a full-wave rectified sine wave, with zeros at the ghost notch 
frequencies and peaks of amplitude 2.0 (6 dB) at the peak frequencies. 

Figure 2 shows the amplitude of a typical ghost filter as a function of frequency. This 
figure shows the ghost filter for the case z = 12m, c = 1500m/s, and with a reflection 
coefficient of -1 .0 at the sea surface. It will be seen that the amplitude decreases to 0 at 
the notch frequencies of 0Hz, 62.5Hz, 125Hz..., and that there are maxima in the 
amplitude at the peak frequencies of 3 1 .25Hz, 93.75Hz... 

Seismic deghosting is an old and long-standing problem. Traditionally, the primary 
goal of deghosting has been to broaden the bandwidth of the data through the notch 
frequencies. Many prior approaches to de-ghosting seismic data have used the 
assumptions of a perfectly flat sea surface, a streamer at a known depth below the sea 
surface, and vanishing pressure at the sea surface. In real seismic data acquisition,' 
however, the sea surface is very often rough. Deterministic algorithms that use an 
explicit mean streamer depth in the deghosting calculations will not work reliably in 
rough sea conditions. Even statistical methods cannot fully correct for the ghost 
reflections in rough sea conditions, owing to the time variant nature of the problem, as 
shown by E.Kragh and R.Laws in "Rough seas and statistical deconvolution", 62 nd 
Annual EAGE Meeting (2001). 

The effects of rough seas on the seismic data have recently been the subject of 
considerable research. In particular, R.Laws and E.Kragh have shown, in "Time-lapse 
seismic and the rough sea wavelet", 70 th Ann. Int. Mtg of Soc. Exploration 
Geophysicists, Extended Abstracts pp 1603-1606 (2000), that errors arising in 
processing data acquired in rough sea conditions are significant for time-lapse seismic 
surveying and for the reliable acquisition of repeatable data for stratigraphic inversion. 

As pointed out in UK patent application No 0015810.5, correction of the rough-sea 
surface ghost requires that the wavefield is completely separated into its up- and down- 
going components. The separation requires in general that the pressure and the vertical 
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component of the particle velocity (or equivalently the vertical pressure gradient) are 
both recorded (see, for example, the equations given by L.Amundsen in Geophysics, 
Vol 58, ppl335-1348 (1993)). 

Another possibility, as suggested by J.O. A.Robertsson and E.Kragh in "Rough sea 
deghosting using a single streamer and a pressure gradient approximation", submitted to 
Geophysics (2001) and further investigated by R#sten et al. in "Rough sea deghosting 
using a vertical particle velocity field approximation", submitted to 63 rd Annual EAGE 
meeting, is to derive an estimate of the vertical component of the particle velocity from 
the pressure acquired at a receiver by using the fact that the pressure vanishes at the sea 
surface. The zero pressure along a contour vertically above the streamer may be 
regarded as a second streamer in a dual streamer configuration. 

Robertsson and Kragh (above) and UK patent application No 0015810.5 suggest using a 
technique from synthetic finite-difference modelling of seismic wave propagation, 
known as the Lax -Wendroff correction, to derive equations for the vertical pressure 
gradient at a streamer in the vicinity of a rough sea surface. The equations express the 
vertical pressure gradient in terms of the pressure gradient along the streamer and the 
time-derivative of the pressure, both of which may readily be obtained from data 
acquired by the streamer. The seismic data may then be de-ghosted using the pressure 
and the approximation for the vertical pressure gradient. In contrast to other proposals 
for seismic deghosting, the method of Robertsson and Kragh is local: a short spatial 
filter (typically of length three) is used to deghost the pressure data acquired at a single 

■ 

receiver. 

A first aspect of the present invention provides a method of determining a vertical 
component of particle velocity from seismic data acquired at a receiver disposed within 
a water column, the method comprising determining the vertical component of particle 
velocity from the pressure acquired at the receiver using a operator accurate for seismic 
data acquired at a depth beneath the surface of the water column of up to at least 0.4 
times the minimum wavelength of seismic energy used to acquire the seismic data. In a 
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preferred embodiment, the operator is accurate for seismic data acquired at substantially 
any depth beneath the surface of the water column. 

The present invention provides an improved method of estimating of the vertical 
particle velocity. Earlier methods of estimating the vertical component of the particle 
velocity have been applicable only for a restricted range for the depth of the receiver 
below the surface of the water column. As discussed below, the prior art method is 
accurate only if the depth of the receiver, below the surface of the water column, is no 
greater than X/3.5, where X is the minimum wavelength of seismic energy used in the 
acquisition process. This severely restricts the choice of receiver arrays with which the 
method can be used. The present invention provides a method of estimating the vertical 
component of the particle velocity that can be used with data acquired at greater depths 
than hitherto. 

The operator may be a spatial and/or wavenumber dependent operator, 

Li preferred embodiment, the method comprises deteimining the vertical component of 
particle velocity from the pressure acquired at the receiver using the following equation 
or an equation derived therefrom: 

v, (a>, x,y,z)= — X F m (k, z)D? (x, y)*p(a> 9 x, y, z) (2) 

P<nz m 

where p is the acquired pressure wavefield, v 2 is the vertical component of particle 
velocity, o is the angular frequency, k co / a) is the wavenumber, a is the P-wave 
velocity, p is the density, D 2 is a differential operator, 

CO 

F m = F m z )= k" lm X A nm fal* 311(1 * is ^ 2D spatial convolution operator. 

n—m 

Equation (2) is an exact expression, and is accurate for any receiver depth. Use of 
equation (2) enables the vertical particle velocity to be obtained from seismic data 
acquired at any receiver depth. 
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In order to reduce the computation required to apply equation (2), it is possible to solve 
an equation derived from equation (2) in which the summation is truncated, rather than 
solving equation (2) exactly. Solving a truncated form of equation (2) reduces the 
computation required, but at the expense of accuracy. A truncated version of equation 
(2) will generally not be accurate for any receiver depth, but will be accurate only for a 
receiver depth less than some maximum value. 

In an alternative preferred embodiment, the method comprises determining the vertical 
component of particle velocity from the pressure acquired at the receiver using the 
following equation or an equation derived therefrom: 

v z (&>, x y z)&F k * p{co, x, z) y where 



min £zj 



b m> a 



n k-(K-i) 



®bnK p 



where an and b m are the backward and forward filter coefficients, Wk are a set of 
positive weights, * is the 1 -D spatial convolution operator and F* = FflcAkx) denote the 
desired discrete horizontal wavenumber response. 

This filter is an approximation of equations (A34) and (A36) given in the appendix 
below. Equations (A34) and (A36) are exact, but are not valid for the case of a rough 
sea surface; this approximation allows them to be applied to the case of a rough sea 
surface. 

In a preferred embodiment, the following filter is used: 

where = cos(kAk x mAx) > b m =0.5g m for g m = 1,2,... ,Af 12 and b 0 =£ 0 . 
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The estimate of the vertical particle velocity provided by the present invention may be 
used to de-ghost the data. Alternatively, the estimate of the particle vertical velocity 
provided by the present invention may be used in other seismic data processing 
applications such as estimation of the source signature. 

A second aspect of the present invention provides a method of processing seismic data 
acquired at a receiver located on a streamer disposed within a water column, the method 
comprising the step of determining a derivative of the pressure in a horizontal direction 
perpendicular to the streamer from a derivative of the pressure in a direction along the 
streamer. This may be done using: 



where the streamer extends along the x-axis and the source is located at (0, ^o). 

In this equation, the receiver co-ordinates in the x-y plane are (xq 9 yo). This enables the 
derivatives of pressure in the cross-line direction to be estimated even for a single 
streamer. In contrast, the prior art methods require at least two laterally separated 
streamers to determine derivatives in the cross-line direction. 

A third aspect of the invention provides an apparatus for determining a vertical 
component of particle velocity from seismic data acquired at a receiver disposed within 
a water column, the apparatus comprising means for determining the vertical component 
of particle velocity from the pressure acquired at the receiver using a operator accurate 
for seismic data acquired at a depth beneath the surface of the water column of up to at 
least 0.4 times the minimum wavelength of seismic energy used to acquire the seismic 
data. 

A fourth aspect of the invention provides an apparatus for processing seismic data 
acquired at a receiver located on a streamer disposed within a water column, the method 
comprising the step of determining a derivative of the pressure in a horizontal direction 



1 a 




dy 



x Q dx 
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perpendicular to the streamer from a derivative of the pressure in a direction along the 
streamer. The apparatus may determine the derivative in the perpendicular direction 
according to: 

ay x 0 ox 

where the streamer extends along the x-axis. 

In a preferred embodiment the apparatus comprises a programmable data processor. 

A fifth aspect of the invention provides a storage medium containing a program for the 
data processor of an apparatus as defined above. 

A further aspect of the present invention provides a method of determining a vertical 
component of particle velocity comprising determining the vertical component of 
particle velocity from the pressure acquired at the receiver using the following equation 
or an equation derived therefrom: 

v, (co> x 9 y 9 z)= — £ F m (k, z)D? (x, y)*p(co, x 9 y, z). 

Preferred features of the invention are defined in the dependent claims. 

Preferred embodiments of the present invention will now be described by way of 
illustrative example with reference to the accompanying drawings in which: 

Figure 1 is a schematic sectional view of a towed marine seismic survey; 

Figure 2 illustrates the effect of ghost reflections on the amplitude spectrum of seismic 
energy emitted by the source in the survey of Figure 1 ; 

Figure 3 is a block flow diagram of a method of the present invention; 
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Figure 4 is a schematic plan view of the seismic surveying arrangement of Figure 1; 

Figures 5(a) to 5(d) illustrate the results of a de-ghosting method of the present 
invention for a streamer at a depth of 6m; 

Figures 6(a) to 6(d) illustrate the results of a de-ghosting method of the present 
invention for a streamer at a depth of 10m; and 



Figure 7 is a schematic block sectional diagram of an apparatus according to the present 
invention. 



Removing the effects of gjiost reflections from pressure data acquired at a seismic 
receiver disposed within a water column is equivalent to separating the acquired 
pressure into its up-going and down-going components. The de-ghosted component of 
the acquired pressure is the up-going component. Various methods for separating 
acquired pressure data into up-going and down-going components have been proposed, 
including the following: 



2 



(1) 



In this equation, P denotes the spatial Fourier transform of the de-ghosted (up-going) 
component of the acquired pressure, P denotes the spatial Fourier transform of the 
acquired pressure p 9 p denotes the density of the water column, go denotes the angular 
frequency, kz is the vertical wave number, and V 2 is the spatial Fourier transform of the 
vertical component of the particle velocity v 2 . 

V zy P and P will in general be functions of the horizontal wave numbers k x and k yy the 
angular frequency, and the depth z of the receiver. The equation (1) above is valid for 
all wave numbers. 
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One prior art method of de-ghosting pressure data acquired at a seismic receiver 
disposed within a water column is to use equation (1) to determine the de-ghosted (up- 
going) pressure component. This method has had the disadvantage that, until now, a 
sufficiently accurate expression for v z has not been available. Although various 
approximate formulations for computing v z have been proposed, these have generally 
only been applicable for a limited range of receiver depth or frequency of seismic 
energy. For example, one prior art approximation for v z is that given in equations (Al) 
and (A2) in the appendix, and this prior art approximation is only accurate if the 
receiver depth is less than approximately X/3.5, where X is the minimum wavelength of 
seismic energy. This approximation also requires that the separation between adjacent 
receivers on a streamer is around 3m or less. Furthermore, this prior method is 
fundamentally limited to frequencies of seismic energy below the first ghost notch 
frequency so that, for example, for a streamer depth of 6m equation (Al) is accurate up 
to only approximately 70Hz. 

The present invention provides more accurate expressions for the vertical particle 
velocity v z , and thus enables more accurate de-ghosting of pressure data using equation 



According to the present invention, the vertical component of the required particle 
velocity is determined using the following equation, or using an equation derived 
therefrom: 



In equation (2), the functions F m (k,z) are functions that are independent of the 
horizontal spatial co-ordinates and are given by: 



(1). 




(2) 



00 




(3) 



WO 03/058281 



PCT/GB03/00049 



11 



where 



(-l)"2 2 "i? 



2n 



(2n)« 



(4) 



In equation (4), Am are scalar quantities, B n is the n Bernoulli number, and n! is the 
factorial of n. 

In equation (2) above D2 is a differential operator, and * denotes the two-dimension 
spatial convolution operator. The components p 9 v 2 are in general functions of a>, x> y 
andz. 

Differentiation generally leads to large amplification of high wavenumbers. Hence, in 
practical situations a bandlimited version of a differential operator is preferably used for 
D 2 . 



The expression for v 2 given in equation (2) is valid for all receiver depths in the water 
column. Use of equation (2) to determine the vertical component of the particle 
velocity enables accurate determination of the vertical component of the particle 
velocity for pressure data acquired at any depth in a water column. This in turn enables 
further processing steps that require knowledge of the vertical component of the particle 
velocity - such as de-ghosting, for example - to be carried out accurately for pressure 
data acquired at any depth in a water column. 



The present invention is not limited to use of the exact equation (2) to obtain the vertical 
component of the vertical particle velocity. The invention also envisages that the 
vertical component of the vertical particle velocity maybe obtained using an equation 
derived from equation (2) such as, for example, an approximate solution to equation (2). 

A special case of equation (2) arises when m is finite, and extends from zero to M. In 
this case, and assuming infinite bandwidth, equation (2) may be re-written as: 
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r a 2 a 
+ 



2 V 



^a* 2 ay 2 J 



(6) 



The functions are given in equation (3) above as infinite series. The functions F m 
may also be written analytically in terms of trigonometric functions, and examples for 
F 0) Fi, F 2 are given in equations (A8), (A9) and (A10) below. 

It will be seen that the x- and y- variables in equation (6) are de-coupled from z. This 
indicates that equation (6) is valid for a rough surface of the water column, as well as 
for a flat surface of the water column. This preferred embodiment of the present 
invention therefore provides a further advantage, since it enables accurate de-ghosting 
of pressure data acquired in rough sea conditions to be carried out, for any receiver 
depth in the water column. 



In the case where M= 2 equation (6) reduces to: 



V « — 

2 



PCDZ 



,2 \ 



dx 2 dy 2 j 



f o2 



,2 Y 



dx 2 dy 



(7) 



The spatial derivatives in equations (6) and (7) may be estimated from the pressure 
acquired at two or more receivers. For example, the pressure derivative along the x- 
axis (which is assumed to coincide with the axis of the streamer) may be determined 
from the pressure recorded at the same time at different receivers along the streamer. 



The preferred method of determining the spatial derivatives will depend on the spatial 
sampling interval - that is, the distance between adjacent sensors on a streamer. Where 
there is a low sampling interval — that is the sensors are closely spaced and so provide 
densely sampled pressure data, the spatial derivatives may be represented by classical 
central finite - difference approximations. In the case of equation (7), for example, the 
second-order spatial derivatives maybe estimated using a three-point finite difference 
operator that estimates the spatial derivative in the x-direction of the pressure at one 
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receiver from the pressure values acquired at that receiver and at the adjacent receiver 
on each side (this is sometimes known as an iterative filter of length 3). 

For a coarse receiver spacing, the horizontal derivatives are most preferably 
implemented as numerically optimised difference operators, or using Fourier 
techniques. In a Fourier technique, the derivatives are obtained by Fourier transforming 
the pressure to the wavenumber domain, multiplying by the negative square of the 
wavenumber (for the second derivative), and inverse Fourier transforming the result. 



Equations (6) and (7) are derived from, and are particular cases of, equation (2), in 
which the variable m is chosen to be finite. Further special cases occur when the 
variables m and n in equation (2) are both chosen to be finite, with m varying from 0 to 
M and with n varying from mXoN (m). In this case, equation (2) reduces to the 
following: 



M 



f a 2 a 2 v 



+ 



{dx* dy') 



(8) 



In equation (8), the functions F m have the same general form as given in equation (3), 
but the summation over n is now carried out from m to N(m) 7 as shown in equation 
(A20) below. 



If we choose M~ 1, N(Q) = 2, and N(l) = 1, equation (8) reduces to the following: 



ptDZ 



( *2 q2 \ 



t j r 



+ 



dx 1 dy 2 J 



(9) 



where the scalar coefficients have the following values: Aoo = 1, Aio = An = -1/3. 

It will be seen that equation (9) has the same general structure as the prior art equation 
(Al) given in the Appendix below. The scalar coefficients Ai 0 and An in equation (9), 
however differ from the scalar coefficients in equation (Al). This has the significant 
effect that equation (9) is valid for greater streamer depths than equation (Al). 
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(Equation (2) above is valid for all streamer depths, but the depth range over which an 
approximate solution to equation (2) will be valid varies from one approximate solution 
to another.) Equation (9) is valid at depths up to approximately A/2. 5, where X is the 
minimum wavelength, whereas equation (Al) requires that the streamer depth is no 
greater than approximately A/3. 5. 

Figure 3 is a block flow diagram illustrating one method of the present invention. 

Initially, at step 9, seismic data are acquired at a plurality of receivers disposed within a 
water column. In general, step 9 will consist of acquiring seismic data using a seismic 
surveying arrangement of the general type shown in Figure 1 in which a plurality of 
seismic receivers SI . . . SN axe disposed on a streamer 2 that is set at a depth h below the 
surface of the water column. The data may be acquired with a single streamer, or they 
may be acquired with two or more streamers that are separated laterally (that is, 
separated along thejy-axis) from one another. 

The invention may alternatively be used to process pre-existing seismic data. Step 9 of 
acquiring seismic data may therefore be replaced by step 10 of retrieving pre-existing 
seismic data from storage. 

At step 1 1 , the vertical component of the particle velocity v z is calculated for one 
receiver (for example the j receiver). Step 1 1 is carried out using equation (2), or 
using any equation derived from equation (2) such as equation (6), (7), (8) or (9). As 
outlined above, the values of the pressure derivatives occurring in these equations are 
estimated using pressure values acquired at two or more receivers. 

Once the vertical component of the particle velocity has been determined for the j th 
receiver, the pressure data acquired at the j* receiver may then be de-ghosted at step 12. 
The de-ghosting step is carried out using equation (1), with the component V 2 being 
obtained by a spatial Fourier transform on the output of step 1 1 . 
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At step 13 the results of step 12 are output, for example by display by an operator, or 
are alternatively stored for subsequent retrieval. 

A receiver Si on the streamer 2 of Figure 1 will acquire data by sampling the incident 
seismic wavefield at discrete time values and so will record the values of the wavefield 
at times to, to + t s , to + 2 t s . . .., where t s is the sampling time. Steps 1 1 to 13 have been 
performed on data acquired at one particular sampling time. At step 14, therefore it is 
determined whether the de-ghosting should be repeated for data acquired by the receiver 
at a different sampling time. If this yields a "yes** determination, a new sampling time 
is selected at step 15, and steps 1 1-13 are repeated for data acquired at the new time. 
This process is repeated until a "no" determination is obtained at step 14. 

At step 16 it is determined whether steps 11,12 and 13 should be repeated for data 
acquired by another receiver on the streamer. If step 16 yields a "yes" determination, a 
new receiver is selected at step 17, and steps 1 1 to 14 and, if necessary step 15, are 
repeated for the new receiver. This loop is repeated until a "no" determination is 
obtained at step 16. 

If it is desired systematically to de-ghost data acquired by every receiver on a streamer, 
the index j may initially be set to 1, step 16 may consist in determining whether the 
index j is equal to the total number of receivers J on the receiver, and step 17 may 
consist of incrementing j by L 

When a "no" determination is obtained at step 16, step 18 then determines whether the 
process should be repeated for data acquired by another streamer. If this step yields a 
"yes" determination, a new streamer is selected at step 19, and steps 1 1-17 are then 
repeated for the receivers on the new selected streamer. 

Where the method of the present invention is applied to data acquired using only a 
single streamer, step 18 maybe omitted and step 19 is not applicable. 
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As noted above, the present invention may be applied to seismic data acquired using 
only a single streamer - that is, to seismic data acquired at receivers all located at the 
same y-co-ordinate value. In this case, it would initially appear difficult to obtain the 
spatial derivatives in the y-direction (perpendicular to the length of the streamer), that 
are required to allow solution of, for example, equations (6), (7), (8) or (9). A further 
aspect of the present invention, therefore, provides a method of estimating the necessary 
derivatives in the y-direction a derivative of the pressure along the streamer. A streamer 
generally contains a large number of pressure sensors disposed at regular intervals along 
the length of the streamer, and the spatial derivative of the pressure along the streamer 
can readily be estimated from the pressure acquired at the sensors. The present 
invention makes it possible to estimate a spatial derivative in a direction perpendicular 
to the streamer from data acquired for a single streamer, from a spatial derivative of the 
pressure along the streamer. 

The method of estimating a spatial derivative in a direction perpendicular to the 
streamer from data acquired for a single streamer may be carried out in connection with 
a determination of the vertical component of the particle velocity according to the 
method described above, but it is not limited to use in connection with determination of 
the vertical component of the particle velocity. 

In a preferred embodiment, a spatial derivative in a direction perpendicular to the 
streamer may be estimated using the following equation: 



The derivation of equation 10 is explained in detail in the Appendix below. 

Figures 5(a) to 5(d) and 6(a) to 6(d) illustrates the results provided by the present 
invention. These results were obtained by modelling synthetic seismic data for a rough- 
sea surface, and de-ghosting the synthetic seismic data using the method of the present 
invention. The synthetic seismic data were computed using the rough-sea modelling 




(10) 
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method developed by R. Laws and E. Kragh (above) and by E. Kragh and R. Laws 
(above). In summary, rough-sea surfaces are generated, and impulse responses beneath 
these surfaces are computed. The impulses are made in sets of three, and a total of 40 
sets were computed with each set being computed from a different realisation of a 4m 
significant wave height sea surface. Data were synthesised for a streamer depth of 6.0m 
(Figures 5(a) to 5(d)) and a streamer depth of 10.0m (Figures 6(a) to 6(d)). 

The synthetic seismic data were then de-ghosted using equation (7) above to compute 
the vertical particle velocity. In the computation using equation (7), the parameters 
were set as follows: M - 1 , N(o) = oo and N(l) = 2 , so giving F 0 = kz cot (kz), 



2 



. 15 ) 



, and F2=0. 



A zero-phase Ricker wavelet with a maximum frequency of 50Hz was convolved with 
the 120 pulses for estimating the v z component. In the estimation of the v z component, 
the well-known length-three central difference operator was used to obtain the second- 
order horizontal spatial derivative in the x-direction. (The calculation was performed 
using a 2-D estimation, so that derivatives in the cross-line direction (y-direction) 
become zero and may be ignored.) 

Once the v z component had been estimated using equation (7), the de-ghosted pressure 
was determined using equation (1) above. The modelled pressure and the estimated v z 
component were down-sampled by a factor of 4 in time before the de-ghosted pressure 
was determined, thereby increasing the temporal sampling distance from 0.25ms to 
1 ,0ms. 



Figure 5(a) shows the raw synthetic seismic data, indicated generally as A, and the 
results of de-ghosting the synthetic data using the Lax-Wendroff technique of 
Robertsson and Kragh (above), indicated as B. 

Figure 5(b) shows the raw synthetic seismic data, indicated generally as A, and the 
results of de-ghosting the synthetic data using equation (7) indicated as C. Figure 5(c) 
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shows the 95% confidence interval in spread of the amplitude spectrum, and Figure 5(d) 
shows the 95% confidence interval in the spread of the phase spectrum. The results are 
shown for the raw synthetic data, indicated as A, the results of de-ghosting the synthetic 
data using the prior art Lax-Wendroff technique, indicated as B, and the results of de- 
ghosting the synthetic data using equation (7) indicated as C. Although the rough-sea 
perturbations are not completely removed from the de-ghosted data obtained by the 
present invention, the correction shows a clear improvement below the first ghost-notch 
frequency (approximately 125Hz). The residual error is of the order of 0.5-1. OdB in 
amplitude and about 5° in phase, and these errors are significantly lower than the errors 
currently obtained in de-ghosting real seismic data. 

Figures 6(a) to 6(d) correspond generally to Figures 5(a) to 5(d) respectively, except 
that they relate to a streamer depth of 10.0m. For a streamer at this depth, the first 
ghost notch frequency is approximately 75Hz. 

A zero-phase band pass filter with cut off frequencies of 10-15-90-100Hz was applied to 
the wavelets shown in Figure 2, whereas a zero-phase band pass filter with cut-off 
frequencies 10-15-65-70 was applied to the wavelets shown in Figure 6(a) to 6(d). This 
has a similar effect to the bandlimited differentiation mentioned above and also below 
in connection with equations (A14) and (A 15). Here the filter is applied to the data 
directly in a pre-processing step, and this has a similar overall effect. 

The de-ghosted seismic data obtained by a method of the present invention may be 
subjected to further processing. For example, the de-ghosted data may be used to 
estimate the signature of the seismic source 1. This can be done, for example, according 
to the method proposed by L. Amundsen et al in Geophysics Vol 60, pp212-222 (1995). 

Alternatively, the de-ghosted data may be further processed to remove or at least 
attenuate events in the data arising from paths that involve reflection and/or scattering at 
the free- surface of the water column. Such events are generally known as free-surface 
related multiples, since they arise from paths that make multiple passes through the 
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water column. In one embodiment, the principal steps required to remove free-surface 
related multiples are as follows. 

Initially the direct arrival and the ghost arrival associated with the direct arrival in the 
incident pressure wavefield acquired at the streamer are isolated and removed from the 
pressure wavefield. This may be done in any suitable way, for example, by "mute" in 
which the direct arrival and its associated ghost arrival is simply separated from the 
other arrivals. Provided that the source-receiver offset is not too great and the water 
depth is not too shallow, the direct arrival in the seismic data, and its associated ghost 
arrival, are generally well separated from other arrivals. 

The vertical component of the particle velocity is then estimated, according to a method 
as described above. The vertical component of the particle velocity is then used to 
separate the remaining incident pressure wavefield (i.e., the pressure wavefield 
remaining after the direct arrival and its associated ghost have been removed) into its 
up-going and downgoing components, for example using equation (1) above to 
determine the up-going component of the pressure wavefield. 

Once the downgoing component of the remaining part of the pressure wavefield has 
been found, it is added to the original direct arrival and its associated ghost arrival. This 
gives the overall downgoing wavefield. (In a typical survey the source will be nearer 
the sea surface than the receiver, so that the direct arrival and its associated ghost 
contain only down-going energy.) Effects arising from free-surface multiple reflections 
may then be removed using the "U divided by D" method proposed by L. Amundsen in 
"Elimination of free-surface related multiples without need of the source wavelet**, 
Geophysics, Vol. 66, pp327-341 (2001). In this method, 'TJ divided by D" refers to the 
up-going component of the wavefield remaining after removal of the direct arrival and 
its associated ghost arrival, divided by the overall down-going component of the 
wavefield. This division is essentially a deconvolution step. 

The embodiments of the invention described above may be carried out on seismic data 
acquired by one or more pressure sensors disposed within a water column, for example 
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on pressure data acquired by a streamer. Alternative embodiments of the invention will 
now be described in which the seismic data include both pressure data and particle 
velocity data. Such data may be acquired by, a seismic receiver, for example a streamer, 
that has both pressure sensors and particle velocity sensors. A streamer provided with 
both pressure sensors and particle velocity sensors can directly measure both the fluid 
pressure p and the particle velocity v. 

The quantities p and v are not independent from one another, but are related by the 
following equations: 



p = -KV.v 



v = -— Vp 
P 



(11) 



(12) 



In equation (1 1), p denotes the time-derivative of p, and K denotes the bulk modulus of 
the water column. 

In one embodiment of the invention, it is assumed that the depth of the streamer is 
known. The depth may be known if, for example, the streamer is provided with depth 
sensors. If the streamer depth is known, it is then possible to compare the output from 
the particle velocity sensors with the output from the fluid pressure sensors. This 
comparison can be made, for example, by calculating the particle velocity from the 
measured values for pressure and from the depth of the streamer according to equation 
(2) or an equation derived therefrom, and comparing the calculated value of the particle 
velocity with directly measured values for the particle velocity. These values should, of 
course, be the same if the particle velocity and the fluid pressure have been measured 
accurately. This embodiment can be considered as calibrating the pressure sensors 
relative to the particle velocity sensors, and allows the accuracy of the sensors to be 
monitored. 



WO 03/058281 



PCT/GB03/00049 



21 

This process for monitoring the accuracy of the sensors can be carried out while 
acquired seismic data is being processed by a method of the present invention. 
Alternatively, it can be carried out on its own, for example to check the sensors at the 
start of a seismic survey. This allows the extra degree of freedom provided by 
measuring both p and v to be used to calibrate the sensors. 



In this embodiment, the depth of the streamer can be estimated, for example from 
knowledge of the arrangement for suspending the streamer from floats. For improved 
accuracy, however, it is preferable that the depth of the receiver is measured directly. 

The invention may also be used to determine the depth of the receiver. To do this the 
vertical component of the particle velocity is determined using a method as described 
above. The depth of the receiver may be estimated from the measured value for the 
vertical component of the particle velocity and the estimated value for the vertical 
component of the particle velocity. 

In an alternative embodiment of the invention, the vertical component of particle 
velocity from the pressure acquired at the receiver is determined using the following 
equation or an equation derived therefrom: 



v z (<y, x 9 z)*F k * p(a>, x, z), 



(13) 



where 



mm 



Fk= b a ^ W " 



(14) 



In equation (14), an and b m denote the backward and forward filter coefficients, Wk are a 
set of positive weights, and F* = F(kAkj) denote the desired discrete horizontal 
wavenumber response. The derivation of equation (14) is explained in the Appendix 
below. 



As noted above, equation (14) is an approximation of equations (A34) and (A36) given 
in the appendix below. Solutions to equation (14) exist for any depths. Equation (14) 
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may be applied to data obtained with a rough sea surface, which is not possible for 
equations (A34) and (A3 6) (although they are exact for the case of a flat sea surface). 
As explained in the Appendix, equation (14) may be implemented to give the following 
zero-phase non-recursive filter: 

~ min{=} , , 2 

Fk= E^ibnSm-^l (15) 

where = cos(kAk x mAx) , b m =0.5g m for g m = l,2,...,M/2 and b 0 =g 0 . 

Vertical particle velocities obtained using equations (13) and (14) or using equations 
(13) and (15) may be used in any of the applications described herein. 

Equations (13), (14) and (15) are 1-D equations, but the corresponding 2-D equations 
may also be used. For the 2-D case, the F k are a function of both k x and ky. . 



The embodiments of the present invention described above have related to the de- 
ghosting the signal recorded at the receiver. The invention can, however, also be used 
for source-side de-ghosting. In principle this could be done by directly measuring the 
fluid pressure at the source, but in practice it is possible to make use of the principle of 
reciprocity to avoid the need to do this. 

The principle of reciprocity is a fundamental principle of wave propagation, and states 
that a signal is unaffected by interchanging the location and character of the sources and 
receivers. For example, if a surveying arrangement with a seismic source at point A and 
a receiver array at point B gives a certain signal at the receiver array, then using a single 
receiver at point A and a source array at point B would lead to the same signal, provided 
that the source array corresponds to the receiver array. (By "corresponds", it is meant 
that the source array contains the same number of sources as the receiver array has 
receivers, and that the sources in the source array are arranged in the same positions 
relative to one another as the receivers in the receiver array). 
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In consequence of the principle of reciprocity, the signal emitted by an array of seismic 
sources can be de-ghosted by making measurements using a receiver array, whether the 

» 

seismic sources are actuated simultaneously or consecutively. The recorded signal is 
analogous to that recorded in a reciprocal arrangement where the source and receiver 
locations are interchanged. The method outlined hereinabove can therefore also be 
applied to the source array that it is desired to de-ghost. The signal produced at the 
receiver array by an array of a plurality of seismic sources is measured, and a de- 
ghosting filter is derived from this measured signal. By the principle of reciprocity, this 
filter can be used to de-ghost the signal emitted by the source array. 

The methods described in this application are fast, and it is possible to use process a 
complete trace using a method of the invention. If desired, however, the methods may 
be applied to a selected part of a trace. 

Figure 7 is a schematic block diagram of an apparatus 20 that is able to perform a 
method according to the present invention. 

The apparatus 20 comprises a programmable data processor 21 with a program memory 
22, for instance in the form of a read only memory (ROM), storing a program for 
controlling the data processor 21 to process seismic data by a method of the invention. 
The apparatus further comprises non-volatile read/write memory 23 for storing, for 
example, any data which must be retained in the absence of a power supply. A 
"working" or "scratch pad" memory for the data processor is provided by a random 
access memory RAM 24. An input device 25 is provided, for instance for receiving 
user commands and data. One or more output devices 26 are provided, for instance, for 
displaying information relating to the progress and result of the processing. The output 
device(s) may be, for example, a printer, a visual display unit, or an output memory. 

The seismic model, the parameters of a seismic surveying arrangement and the prior 
AVO information may be supplied via the input device 25 or may optionally be 
provided by a machine-readable data store 27. 
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The results of the processing may be output via the output device 26 or may be stored. 

The program for operating the system and for performing the method described 
hereinbefore is stored in the program memory 22, which may be embodied as a 
semiconductor memory, for instance of the well known ROM type. However, the 
program may well be stored in any other suitable storage medium, such as a magnetic 
data carrier 22a (such as a "floppy disk") or a CD-ROM. 



Appendix 



Roberts son and Kragh (above) have shown that the vertical component of the particle 
velocity, v x , can be estimated from the raw measurements of the pressure field, p . In 
the typical seismic frequency band one possible approximation given is 



v z (cD,x,y,z)= - 



pcoz 



A' 00 + 4(feJ + A[ x z 



r *2 



d* d 
+ 



2 \ 



dx 2 dy 



J _ 



(Al) 



with scalars A' Q0 = 1 ; A' l0 =-1/2 ; ^ = -1/2, 



(A2) 



In equation (Al) a> is the (angular) frequency, k = a)/a is the wavenumber, a is the P- 
wave velocity, p is the density, and z = z(x 9 y) is the streamer depth. According to 
Robertsson and Kragh, approximation requires that the streamer is towed no deeper than 
approximately A/3.5 , where X is the minimum wavelength, and that the receiver 
spatial sampling interval is around 3m or less. This method is fundamentally limited to 
frequencies below the first ghost notch. For instance, for a streamer depth of 6 m, 
equation (Al) is accurate up to approximately 70 Hz. 

For densely sampled pressure measurements one attractive feature of method (Al) is 
that the filters for estimating the vertical particle velocity are local and short. By 
approximating the second-order horizontal spatial derivatives by three-point central 
difference-operators, the spatial filter has cross structure with five coefficients. 
Processing single streamer pressure data such as those provided by the new Q-Marine 
streamers under the 2.5D earth model, the filter has three coefficients only. Thus, to 
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estimate the vertical particle velocity at a particular receiver location, only the pressure 
field at this and the two neighbouring receivers are required. Provided the depth of the 
(near-horizontal) streamers below the sea surface are known at all points along the 
streamers, equation (Al) gives a methodology of estimating the vertical particle velocity 
field in rough-sea environments . Let P and V x denote the spatial Fourier transforms of 
p and v 2 , respectively. The estimate given in equation (Al) then can be combined 

with the raw pressure measurements to deghost the data (see e.g. L. Amundsen in 
Geophysics Vol. 56, ppl027-1039 (1993) or UK patent application No. 9906456,0) 



P(p,k x ,k y ,z)=^ P(p,k x ,k y ,zy^V t (p,k xi k y9 z) 



K 



(A3) 



where p{a>,k x ,k y9 z) denotes the deghosted (up-going) pressure, k z = V& 2 -k 1 is the 

2 2 2 

vertical wavenumber, k -k x + k y , and k x and k y are the horizontal wavenumbers. 

Equation (3) is valid for all wavenumbers. Alternatively, approximations to equation 
(A3) as those suggested by A.Osen et al. in "Towards optimal spatial filters for 
demultiple and wavefield splitting of ocean bottom seismic data", submitted to 
Geophysics (2002). 

In the present application the accuracy of the vertical particle velocity field 
approximation given in equation (Al) is extended by introducing optimised filters. 
L.Amundsen et al. derived, in Geophysics Vol. 60, pp2 12-222 (1995) an exact 
frequency-wavenumber domain relationship between the vertical component of the 
particle velocity and the pressure field. For single streamer data the relationship is: 



V,L,k„k„z)-±- eX T^T eX T^ kk*,.*,.*) (A4) 



'exp(- ik z z)+ exp^z)^ 



pco 



exp(- ik z z)- exp(/^z) j 



Equation (A4) is valid for all streamer depths z in the water column when the receiver 
cable is horizontal and the air/water surface is flat with vanishing pressure. Essentially, 
the equation shows that to find the vertical component of the particle velocity from the 
pressure field, first the pressure has to be receiver-deghosted, then the receiver ghost 
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operator for the vertical component of the particle velocity is applied. (Taking into 
account that the reflection coefficient of the free surface is -1, the ghost operator for the 
pressure is 1 - exp{2ik t z) since the pressure is the sum of upgoing and downgoing 

waves, whereas the ghost operator for the vertical particle velocity is 1 + exp(2ik 2 z) 

since the vertical particle velocity is proportional to the difference between upgoing and 
downgoing waves.) Note that when the source depth is less than the receiver depth the 
incident pressure wavefield (the sum of the direct pressure wavefield and its ghost) does 
not contain a receiver ghost, but only a source ghost. Since equation (A4) relies on 
receiver deghosting, it cannot correctly treat the incident wavefield. This is of minor 
importance in receiver deghosting of streamer data. When the source depth is greater 
than the receiver depth, however, the incident pressure wavefield contains a receiver 
ghost as do the reflected part of the pressure field. In this case, equation (A4) is valid 
for the full pressure wavefield. 

Equation (A4) can be rewritten as 



GO 



with horizontal wavenumber ( k ) independent functions 



F m = F m Qc,z)= k-^A m (kzr (A6) 



where 



(-l)"2 2 -5 3 



m ~ (in) 



(A7) 



are scalars, B H is the n:th Bernoulli number, and n\ is the factorial of n . All functions 
F m can be given in compact form. For instance, 



F 0 = Azcot(Az)= —ikz 



f exp(- ikz) + exp(/fe) 
^exp(- ikz) - exp(ifcz) y 



(A8) 
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A- Z 



2 4* 3 



£z-2 



d(kz) ) d(jcz) 



d 



F 0 . (A10) 



The fixnction F 0 ensures that vertically propagating waves (k x —k y = 0) are correctly 
deghosted. 

Taking an inverse Fourier transform of equation (A5) over horizontal wavenurnbers 
gives 



v r (o>, x,y, z) = — i- £ F m (fc, z)D a " (x, j/)*p(fi>, x, z), (All) 

m-0 

where Z> 2 is the inverse Fourier transform of tc 2 , and * is the 2D spatial convolution 
operator. For infinite bandwidth, 

D 2 (x,y)=D 2 (x)+D 2 (y), * (A12) 

where 



"£r ; D,(y)=-^, (A13) 

are second-order spatial differentiation operators. Differentiation however involves 
large amplification of high wavenurnbers. Hence, in practical situations a bandlimited 
version of differentiation is preferably used. Let d 2 denote a bandlimited version of 

D 2 . Two possible bandlimited differentiation operators are 



d 2 (*)= — + Ax)- 2<?(x)+ S(x - Ax)], (A14) 
Ax 



where Ax is the spatial sampling interval, or 
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d 2 (x,co)= - \dk x oos(k x x)v(a> 9 k x %k x 2 ) (A15) 

n i 

where W is a weighting function that properly bandlimits the differentiation operator. 
Note that W can be changed per frequency component. In equation (A15) K is the 
maximum wavenumber. Two possible choices are K = co/c 9 where c is the velocity of 
water, or K = n I Ax , the Nyquist wavenumber. 

Note that equation (Al 1) is valid for all receiver depths z in the water column. Also, 
observe that the scheme is iterative: 

D?*p = D?*(D?- l *p) ; £> 2 °=1. (A16) 

This recursive property allows equation (Al 1) to be implemented with numerical 
efficiency. For dense spatial sampling the differentiators D 2 (x) and D 2 (y) can be 
approximated by filters of length three [see equation (A 14)]. A properly designed, 
bandlimited, second-order differentiator d 2 (x) is imperative for obtaining a stable 
iterative scheme. 

Special case I: Finite m 

A special case of equation (Al 1) arises when m is finite, say M . Then, for infinite 
bandwidth, 



m 



v M (p 9 x t y 9 z)* — Z F «fc z i"ZT + T7 Pfa> x >y> z )> ( A1? ) 



where F m is given in equation (A6) as an infinite series. However, the function F m can 
also be given analytically in terms of trigonometric functions. Examples are F 0 , F l9 
and F 2 given in equations (A8), A(9), and (A10), respectively. Observe that the 
assumption of flat sea surface used to derive equation (4) can now be relaxed. In 
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equation (A17) the x and y variables are decoupled from z , implying that it is valid 
also for a rough sea surface. 

Equation (A17) can be numerically implemented in different ways. The preferred 
implementation depends on the spatial sampling interval for the pressure field. For 
densely sampled data one would typically represent the spatial derivatives by classical 
central finite-difference approximations. If M is small, the filter is then local and 
compact. As an example, choosing M = 2 gives for infinite bandwidth, 

p((D 9 x 9 y,z). (A18) 

Approximating the second-order spatial derivatives by a three-point finite difference 
operator the iterative filter is of length three. Since the functions F Q , F l and F 2 are 

valid for all receiver depths, equation (A18) is more general than equation (Al), and is 
preferred for rough-sea deghosting. 

For a coarser receiver interval the horizontal derivatives are most conveniently 
implemented as numerically optimised difference operators or by use of Fourier 
techniques. 



v z (p) 9 x,y 9 z)tt 



pcoz 



a 2 a 
+ 



2 > 



dx 2 dy 2 ) 



r ^2 



d* d 



2 \ 



v 



dx 2 dy 2 ) 



Special case II: Finite m, n 

Limiting the sums over m in equation (11) and n in equation (A6) from infinity to M 
and N(m) , respectively, gives for infinite bandwidth the approximation 



j* M J Q 2 Q 2 V 

v r (co, x, y,z)* 2 F m (k, z) — + — p(o>, x y y, z), 



pox—* 



dx' dy* ) 



(A19) 



with 



m 



(A20) 



n*=m 
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Choosing M = 1 , N(o) = 2 , and N(l) = 1 , gives 



(A21) 



(A22) 



and 



V s (fl>,^y,z)w - 



f *2 



3 2 a 



2 \ 



^&c 2 0y 2 J 



p(tf>,*,.y,z) 9 



(A23) 



with scalars 



4>o=l ; 4 0 =-l/3 ; ^=-1/3. 



(24) 



Equation (A23) has the same structure as equation (Al) proposed by Robertsson and 
Kragh (above). The scalars A i0 and A n , however, differ from the scalars A[ Q and A[ x 

in equation (Al). In many respects the limitations of equation (A23) are similar to 
those of the formula of Robertsson and Kragh (above). However, equations (Al 1), 
(A17), and (A19) are general and one significant difference is that they are valid for all 
streamer depths, whereas the theory derived by Robertsson and Kragh (above) requires 
that the streamer is towed no deeper than approximately A/3.5 , where X is the 
minimum wavelength. 



If we choose M = 1 , N(o) = oo and N{\) = 2 , this results in F 0 = kz cot(fe) and 



r 2 N 
1 + — k 2 z 2 

15 J 



V 



Data processing of single streamer data 



Equations (Al 1), (A17) and (A19) above include derivatives in the ^/-direction or cross- 
line direction (the streamer(s) is/are assumed to extend along the x-direction). In a 
seismic surveying arrangement that includes only a single streamer it is not clear how 
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the cross-line derivatives in the y -direction that occur in the above equations may be 

determined from the seismic data acquired by the receivers on the streamer. It will be 
noted, however, that only even powers (second-order, fourth-order, sixth-order, etc.) 
spatial derivatives occur. Finding an expression for the second-order spatial derivative 
is therefore sufficient. 



Consider the single streamer geometry shown in Figure 4, which is a schematic plan 
view of the seismic surveying arrangement of Figure 1. The seismic source 1 is located 
at coordinates (0,j>o) in the x-y plane. The streamer extends along the *-axis, and has a 
y-coordinate yo. One of the receivers Si is located at coordinates (xojo), and it assumed 
that the source 1 and the receiver S\ are located in the same z-plane in an acoustic layer. 
It will be assume that the model only varies with depth, and does not vary with x oxy. 
By using a finite-difference approximation to a second-order derivative, it is possible to 
write the following equation for the pressure p : 

p( x o>yo)=-rr(p(xo>yo + &y)- 2 p( x o>y Q )+/>(w<> -A^))+o^y 2 ), (A25) 



dy^^ Ay 

where o[Ay 2 ) expresses the leading error term in the finite-difference approximation. 

However, because of radial symmetry the recorded pressure is the same on either side of 
the receiver so that equation (A25) becomes 



p(? 0 ,y 0 )= — t(p(w<> + Ay)-/>(x 0 ,jy 0 ))+O^y 2 ). (A26) 



dy" ^ Ay 
Radial symmetry also yields 



p(x 0 , y 0 + Ay)= p\x 0 -Jl + (Ay/x o y, y 0 \ (A27) 



so that equation (A26) becomes: 



P(x 0 ,y 0 )= —2 [p[x 0 yjl + (Ay I x 0 J , y 0 J - p(x 0 , y 0 )J + o(Ay 2 ) (A28) 



dy^ Ay 
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By using a first-order one-sided finite-difference approximation in the x -direction it is 
possible to write equation (A28) as: 



dy 



(Wo )= " — 2 l \j£p(* 0 ,y 0 )+ o[Vl + (Ay/xJ- ljj + 0(&y 2 ) 



(A29) 



Taking a Taylor series expansion of the square root: 
<Jl + (fr/xJ = l + ^(Ay/xJ+o(Ay 4 ), 



(A30) 



equation (A29) may be re-written as , 



2x, 



dy 



p(x 0 ,y 0 )= 



^(Ay/xJ + O^) 



J 



Ay 



f d_ 



p(x 0 >y 0 )+o(Ay 2 j\ + o(A.y 2 ): (A31) 



Further simplification of equation (31) yields 



p( x <> > y* )= — p( x o > y 0 )+ o{Ay 2 ). 



dy 



x 0 dx 



(A32) 



Finally, we let Ay -> 0 so that: 



9 p( x o>y 0 )=—4-p(x 0 ,y 0 ). 



dy 



x 0 dx 



(A33) 



For pressure data over a horizontally layered medium there is no problem with the 

singularity in equation (33) since — _p(*o>.Vo) 6 oes to zero faster than the singularity 

dx 



1 / x 0 goes to infinity. 
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Equation (A3 3) is consistent with what is valid for cylindrical symmetry (data over a 
horizontally layered medium). An alternative approach would therefore be to multiply 

the data with , where t is time, and use the 2D version of the theory derived in this 
application as well as in UK patent application No. 0015810.5. However, equation 
(A33) can potentially give more accurate results. 

It should be noted the above approach may not be reliable for a medium that is not 
simply plane layered since there can be a problem because p is not symmetric with 

respect to x. 



Complex frequency 



In the implementation of the filters we use complex frequencies to avoid the 
singularities in the functions F m as given e.g. in equations (8), (9), and (10) for 

M = 0,1,2 , respectively. The use of a complex-frequency is described by, for example, 

S.Mallick and L.N.Frazer in Geophysics Vol. 52, ppl355-1364 (1987) and by 
L.Amundsen and B. Ursin in Geophysics Vol. 56, ppl027-1039 (1991). 

On V 7 -estimation using numerically optimised spatial filters 



1. Introduction 

In this section, we introduce numerically optimised spatial filters for the estimation of 
the vertical particle velocity field. Amundsen et al. (1995, above) derived an exact 
frequency-wavenumber domain relationship between the vertical component of the 
particle velocity and the pressure field. For single streamer data the relationship in the 
2D case is 



v,(<a,A: x ,z) = - 



exp(—ik x z) + exp(ik z z) 



pco (_exp(— ik z z) - exp(ik z z) 



p(a),k x ,z) 



(A34) 



Equation (A34) may be re-written as: 
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v z (to, k x ,z) = F(co, k x z) ■ pip, k x , z) (A35) 



where the v 2 -estimation filter is given by: 



exp(-ik 2 z) + exp(zfc t z) 
pco |_exp(-z£ x z) - exp(ik z z) 



In the frequency-space domain, equation (A34) is written symbolically as 



(A36) 



K z = F*p {All) 



where * denotes spatial convolution. 



2. Optimisation problem formulation 



The discrete horizontal wavenumber response of a spatial digital filter can be written as: 



M 12 

y^ j b m cxp(-ikAk x mAx) 

F k = F(kAk x ) = = m (A38) 

^a n exp(-ikAk x nAx) kn °' fX 

where a n> b m are the backward and forward filter coefficients with N,M the respective 

filter orders. Without loss of generality, M and N are assumed to be even numbers 
and a 0 si, Furthermore, Ax is the spatial sampling interval and Ak x = 7u/\Ax(K -1)] 

is the distance between the discrete horizontal wavenumbers. If F k = F(kAk x ) denotes 

the desired discrete horizontal wavenumber response, then the design of the spatial 
filters can be stated as a general weighted non-linear least-squares optimisation 
problem: 
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Ft = 



mm 



K-i 



®^b 

km m Jp* 



(A39) 



where W k is a set of positive weights. In the following, the optimisation problem is 

formulated for propagating waves only, for which the ideal decomposition filter is 
purely real (i.e., zero-phase). Let k denote the index corresponding to a discrete 
horizontal wavenumber close to the critical wavenumber. 



3. Zero-phase FIR filters 



For a zero-phase non-recursive or finite impulse response (FIR) filter, the forward filter 
coefficients are symmetric and the backward filter order is zero. This reduces the non- 
linear least-squares optimisation problem given in equation (A3 9) to the following 
linear-least squares optimisation problem: 



mm 



K-l 



2 



Z w *Ks»- F *\ ( A4 °) 

where = oos(kAk x mAx) . The g m 's are related to the b m 's by b m - 0.5g m for 
g m =l,2,...,M/2 and by b Q =g 0 . 



4. Quadratic Programming 



The linear least-squares optimisation problem shown in equation (A40) is 
unconstrained, and unique and optimal solutions are given by the pseudo-inverses. 
Nevertheless, it can be important that the magnitude responses of the optimal FIR filter 
have zero-error at particular horizontal wavenumbers, or that they have magnitude 
responses smaller than a given function for a certain horizontal wavenumber range. The 
combination of a linear least-squares optimisation problem with linear equality and/or 
inequality constraints is called Quadratic Programming (QP) and can be solved 
efficiently with many different algorithms. 
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For example, a single equality constraint can be included at zero horizontal 
wavenumber to ensure that propagating waves are estimated correctly by the optimal 
FIR filter. In addition, linear inequality constraints can be included in the evanescent 
horizontal wavenumber range to force the spatial filter's magnitude response of the 
ideal decomposition filter. This prevents boosting of evanescent waves by the 
numerically optimised spatial filters. 
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CLAIMS: 

1 . A method of determining a vertical component of particle velocity from seismic 
data acquired at a receiver disposed within a water column, the method comprising 
determining the vertical component of particle velocity from the pressure acquired at the 
receiver using a operator accurate for seismic data acquired at a depth beneath the 
surface of the water column of up to at least 0.4 times the minimum wavelength of 
seismic energy used to acquire the seismic data. 

2. A method as claimed in claim 1 wherein the operator is accurate for seismic data 
acquired at substantially any depth of the receiver beneath the surface of the water 
column. 

3. A method of determining a vertical component of particle velocity as claimed in 
claim 1 and comprising determining the vertical component of particle velocity from the 
pressure acquired at the receiver using the following equation or an equation derived 
therefrom: 

■ 

v z fa, x, y y z)= — £ F m (k 9 z)D? (x, y)*p(cD, x 9 y, z), 

where a> is the angular frequency, k( = col a) is the wavenumber, a is the P-wave 
velocity, p is the density, D 2 is a differential operator, 

co 

F m ~ F m (k, z)= k~ 2m ^] (kzj n and * is the 2D spatial convolution operator. 

n«=m 

4. A method of determining a vertical component of particle velocity as claimed in 
claim 2 and comprising determining the vertical component of particle velocity from the 
pressure acquired at the receiver using the following equation: 



(<i>, x, y, z)= — J] F m (k, z)D 2 M (x, y)*p(p, x, y, z), 
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5. A method as claimed in claim 3 or 4 wherein D 2 is a bandwidth limited 
differential operator. 

6. A method as claimed in claim 3 wherein the summation over m is carried out for 
m = 0 to m = M, 

7. A method as claimed in claim 6 wherein the step of determining the vertical 
particle velocity comprises determining the vertical component of particle velocity from 
the pressure acquired at the receiver according to 

p(p) 9 x 9 y,z). 



r *2 



pcaz 



2 \ 



+ 



{dx* dy 



r *2 



J 



d* a 
+ 



2 Y 



{dx 2 dy 2 ) 



8. A method as claimed in claim 6 wherein the step of determining the vertical 
particle velocity comprises determining the vertical component of particle velocity from 
the pressure acquired at the receiver according to: 



i M 4 d 2 d 2 V 



dx L dy L ) 



where F m & k 

HI 



-2 m 



AM" 



9. A method as claimed in claim 6 wherein the step of determining the vertical 
particle velocity comprises determining the vertical component of particle velocity from 
the pressure acquired at the receiver according to: 

v t (p,x,y t z)* — 

pcaz 

where -^oo = 1, and A\§ —A\\ = -1/3. 



+ 



dx l by 1 ) 



p(p,x y y y z)> 
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10. A method as claimed in claim 3 and comprising determining the vertical 
component of particle velocity from the pressure acquired at the receiver using the 
following equation or an equation derived therefrom: 

v 2 (co, x, y, z) « F k p(o>, x, y y z) , where 



~ min £zj 



where an and b m are the backward and forward filter coefficients, are a set of 
positive weights, and F k = F(kAkjJ denote the desired discrete horizontal wavenumber 
response. 

11. A method as claimed in claim 1 0 and comprising determining the vertical 
component of particle velocity from the pressure acquired at the receiver using th^e 
following equation: 



S m k=0 

where = cos(kAk x mAx) , b m = 0.5g m for g m = 1,2,..., Ml 2 and b 0 = g 0 



12. A method as claimed in any preceding claim wherein the data has been acquired 
at a receiver disposed on a steamer, the method further comprising the step of 
determining a derivative of the pressure in a horizontal direction perpendicular to the 
streamer from a derivative of the pressure in a horizontal direction along the streamer. 

13. A method of processing seismic data acquired at a receiver located on a streamer 
disposed within a water column, the method comprising the step of determining a 
derivative of the pressure in a horizontal direction perpendicular to the streamer from a 
derivative of the pressure in a horizontal direction along the streamer. 
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14. A method as claimed in claims 12 or 13 wherein the step of determining a 
derivative of the pressure in a horizontal direction perpendicular to the streamer 
comprises using: 



where the streamer extends along the x-axis. 

15. A method as claimed in any of claims 1 to 14 and comprising the further step of 
processing the seismic data using the vertical component of the particle velocity. 

16. A method as claimed in claim 1 5 and comprising processing the seismic data so 
as to determine at least one of an up-going component and a down-going component of 
the seismic data. 

17. A method as claimed in claim 1 5 and comprising processing the seismic data 
using the vertical component of the particle velocity thereby to attenuate effects in the 
processed data of seismic energy reflected and/or scattered at the free-surface of the 
water column. 

18. A method as claimed in claim 1 5 and comprising 

a) removing the direct arrival and its associated ghost arrival from the acquired 
seismic data; 

b) separating the seismic data remaining after removal of the direct arrival and its 
associated ghost arrival into its up-going and down-going components; 

c) summing the down-going component obtained in step(b) and the direct arrival 
and its associated ghost arrival thereby to obtain the overall down-going component; 
and 

d) dividing the up-going component obtained in step (b) by the overall down-going 
component obtained in step (c) thereby to attenuate the effect of multiple reflections in 
the seismic data. 



dy 



r-p(wo)= 



x Q dx 
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19. A method as claimed in claim 15 and comprising processing the seismic data 
using the vertical component of the particle velocity thereby to obtain information about 
the signature of the source of seismic energy. 

20. . A method as claimed in any of claims 1 to 14 wherein the seismic data 
comprises pressure and particle velocity data, and wherein the method further comprises 
comparing the vertical component of the particle velocity determined from the acquired 
pressure with the measured values of the particle velocity. 

21. A method as claimed in any of claims 1 to 14 wherein the seismic data 
comprises pressure and particle velocity data, and wherein the method further comprises 
the step of determining the depth of the seismic receiver within the water column from 
the measured vertical component of the particle velocity and the vertical component of 
the particle velocity determined from the pressure. 

22. A method as claimed in claim 7, 8 or 9, or in any of claims 10 to 21 when 
dependent directly or indirectly from one of claims 7,8 or 9, and comprising 
determining at least one horizontal derivative of the pressure wavefield according to a 
method including the step of Fourier transforming the pressure wavefield to the 
wavenumber domain. 

23. A method of acquiring marine seismic data comprising the steps of: 
actuating an array of one or more seismic sources to emit seismic energy; 
acquiring seismic data at one or more receivers disposed within a water column, 

the seismic data including at least pressure data; and 

processing the acquired pressure data according to a method as defined in any of 
claims 1 to 22. 

24. An apparatus for determining a vertical component of particle velocity from 
seismic data acquired at a receiver disposed within a water column, the apparatus 
comprising means for determining the vertical component of particle velocity from the 
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pressure acquired at the receiver using a operator accurate for seismic data acquired at a 
depth beneath the surface of the water column of up to at least 0.4 times the minimum 
wavelength of seismic energy used to acquire the seismic data. 

25. An apparatus as claimed in claim 24 and adapted to determining the vertical 
component of particle velocity using a operator accurate for seismic data acquired at 
substantially any depth of the receiver beneath the surface of the water column 

26. An apparatus as claimed in claim 24 or 25 and comprising means for 
determining the vertical component of particle velocity from the pressure acquired at the 
receiver using the following equation or an equation derived therefrom: 

v 2 (a>, x 9 y 9 z)= — J] F m (k 9 z)D 2 m (x 9 y)*p(a> 9 x 9 y 9 z) 9 

where co is the angular frequency, A: is the wavenumber, a is the P-wave velocity, p is 

QO 

the density, D 2 is a differential operator, F m = F m (k 9 z)- k~ 2m ^ A nm (kzj n and * is the 
2D spatial convolution operator. 

27. An apparatus as claimed in claim 26 and adapted to determine the vertical 
component of the particle velocity using a bandwidth limited differential operator for 
D 2 . 

28. An apparatus as claimed in claim 26 or 27 and adapted to carry out the 
summation over m from m = 0 to m = M. 



29. An apparatus as claimed in claim 28 and adapted to determine the vertical 
component of particle velocity from the pressure acquired at the receiver according to 



v s (a) 9 x 9 y 9 z)x 



pcoz 



+ 



I a* dy'J 



{dx 2 dy 



J 



p(co 9 x 9 y,z). 
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30. An apparatus as claimed in claim 29 and adapted to determine the vertical 
component of particle velocity from the pressure acquired at the receiver according to: 



> a fa, x, y,z)* — — ]T F m fa z \^r + 4-r I p(p> x > y> z )> 

poz^ \dx 2 dy J 



where F„ m k^^A^^T ■ 

n—m 



31. An apparatus as claimed in claim 28 and adapted to determine the vertical 
component of particle velocity from the pressure acquired at the receiver according to: 



v,(fi>,x,j/,z)« 



poyz 



+ 



v 



dx 2 dy 2 ) 



pQo 9 x 9 y 9 z) 9 



where Aq 0 = 1, and^io = = -1/3. 



32. An apparatus as claimed in any of claims 25 to 31 for processing data acquired 
at a receiver disposed on a steamer and further adapted to determine a derivative of the 
pressure in a horizontal direction perpendicular to the streamer from a derivative of the 
pressure in a direction along the streamer. 

33. An apparatus for processing seismic data acquired at a receiver located on a 
streamer disposed within a water column, the method comprising the step of 
determining a derivative of the pressure in a horizontal direction perpendicular to the 
streamer from a derivative of the pressure in a direction along the streamer. 

34. An apparatus as claimed in claims 32 or 33 and adapted to determine a 
derivative of the pressure in the direction perpendicular to the streamer according to: 



oy x Q ox 
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■ 

where the streamer extends along the x-axis. 

35. An apparatus as claimed in any of claims 25 to 34 and adapted to further process 
the seismic data using the vertical component of the particle velocity. 

36. An apparatus as claimed in claim 35 and adapted to process the seismic data 
using the vertical component of the particle velocity thereby to reduce effects in the 
processed data of seismic energy reflected and/or scattered at the surface of the water 
column. 

37. An apparatus as claimed in claim 35 and adapted to process the seismic data 
using the vertical component of the particle velocity thereby to obtain information about 
the signature of the source of seismic energy. 

38. An apparatus as claimed in any of claims 25 to 37 and comprising a 
programmable data processor. 

39. A storage medium containing a program for the data processor of an apparatus 
as defined in claim 38. 
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Method I - Nominal streamer depth of 6.0 m 
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FIG.5(a) to 5(d): Rough-sea deghosting results for a nominal 
streamer depths of 6.0 m.Method I is the original method of Robertsson 

« 

and Kragh (2001) whereas Method II is the new method. 
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FIG.6(a) to 6(d): Rough-sea deghosting results for a nominal 
streamer depths of 10.0 m.Method I is the original method of Robertsson 
and Rragh (200 1) whereas Method EE is the new method. 
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